"""
统计单位面积下有丝分类数量
"""
import os
import numpy as np
import json

xls_name = "sm" + ".xls"
geo_json_folder = "/media/hsmy/16TB/20241110GIST_南方医_中六_核分裂/sm/241130_geojson/"
output_folder = "/media/hsmy/16TB/20241110GIST_南方医_中六_核分裂/sm/mc"
os.makedirs(output_folder, exist_ok=True)
xls_list = []
radius = 5047
# radius = 6300
for geo in os.listdir(geo_json_folder):
    geo_path = os.path.join(geo_json_folder, geo)
    json_data = json.load(open(geo_path))
    features = json_data['features']
    points = []
    for feat in features:
        coordinates = feat['geometry']['coordinates']
        point = coordinates[0][0]
        points.append(point)
    if len(points) == 0:
        continue
    points = np.array(points)
    distances = np.sqrt(np.sum((points[:, np.newaxis] - points) ** 2, axis=2))
    counts = np.sum(distances <= radius, axis=1)
    max_index = np.argmax(counts)

    xls_list.append({
        "file": geo.replace(".geojson", ""),
        "xy": points[max_index],
        "xy_count": counts[max_index],
        "sum": len(points)
    })

    mc_max_x, mc_max_y = points[max_index]
    temp = {"type": "feature", "geometry": {"type": "Polygon", "isEllipse": True, "coordinates": []},
            "properties": {"object_type": "annotation"}}
    y_min = int(mc_max_y - radius)
    y_max = int(mc_max_y + radius)
    x_min = int(mc_max_x - radius)
    x_max = int(mc_max_x + radius)

    a = [x_min, y_min]
    b = [x_min, y_max]
    c = [x_max, y_max]
    d = [x_max, y_min]
    temp["geometry"]["coordinates"] = [[a, b, c, d, a]]
    temp["properties"]["color"] = [255, 0, 0]
    json_data['features'].append(temp)

    output_path = os.path.join(output_folder, geo)
    with open(output_path, 'w', encoding='utf-8') as fw:
        json.dump(json_data, fw, indent=4, ensure_ascii=False)

with open(xls_name, 'w+', encoding='utf-8') as fw:
    head = "file\txy\txy_count\tsum\n"
    fw.write(head)
    for i in xls_list:
        fw.write(f"{i['file']}\t{i['xy']}\t{i['xy_count']}\t{i['sum']}\n")